

%% "true" DGP point estimates
load('../true_DGP/tDGP.mat','DGP')
true_beta=DGP(1:(2*N));
clear DGP


%% reset rng seed
rng('default')


%%
NSP=1e+3; % number of simulated sample paths
T_chn=4;
epsilon=mvnrnd(zeros(N,1),Sigma,(T-(54-T_chn-1))*NSP)'; % growth shocks
epsilon=reshape(epsilon,N,T-(54-T_chn-1),NSP);


%% small shock (6.5%) to Mongolia, Taiwan, India, Japan, South Korea, Philippines, (2004-T_chn):2004
P_chn=zeros(size(P(:,:,54-T_chn:T-1)));
parfor nsp=1:NSP
    ep=epsilon(:,:,nsp);
    y_chn=y;
    y_chn([131:132,134:136,148],54-T_chn:54)=0.065+ep([131:132,134:136,148],1:T_chn+1);
    PP_chn=zeros(size(P(:,:,54-T_chn:T-1))); %#ok<*PFBNS>
    BB_chn=zeros(2*N,T-54+T_chn);
    % 2001
    nm=find(M(:,54-T_chn)==1);
    DD=[diag(1-D(nm,54-T_chn)),diag(D(nm,54-T_chn))];
    PP_chn(:,:,1)=P(:,:,54-T_chn);
    PP_chn([nm;N+nm],[nm;N+nm],1)=PP_chn([nm;N+nm],[nm;N+nm],1)+DD'*(Sigma(nm,nm)\DD);
    Pd_chn=sqrt(diag(PP_chn(:,:,1)\eye(2*N)));
    BB_chn(:,1)=B(:,54-T_chn);
    BB_chn([nm;N+nm],1)=BB_chn([nm;N+nm],1)+PP_chn([nm;N+nm],[nm;N+nm],1)\(DD'*(Sigma(nm,nm)\(y_chn(nm,54-T_chn)-DD*BB_chn([nm;N+nm],1))));
    U0_chn=exp(alpha(130,1)+theta(1)*(Pd_chn(130,1)*sgi_nodes(:,1)'+BB_chn(130,1)+...
        ep_sd(130,1)*sgi_nodes(:,2)'));
    U0_chn=(U0_chn./(1+U0_chn))*sgi_weights;
    U1_chn=exp(alpha(130,1)+theta(2)*(Pd_chn(N+130,1)*sgi_nodes(:,1)'+BB_chn(N+130,1)+...
            ep_sd(130,1)*sgi_nodes(:,2)')-f(130,1));
    U1_chn=(U1_chn./(1+U1_chn))*sgi_weights;
    DD_chn=D;
    DD_chn(130,54-T_chn+1)=M(130,54-T_chn+1).*(U1_chn>=U0_chn);
    for t=1:T-54+T_chn-1
        if DD_chn(130,54-T_chn+t)~=D(130,54-T_chn+t)
            y_chn(130,54-T_chn+t)=true_beta(130)*(1-DD_chn(130,54-T_chn+t))+true_beta(N+130)*DD_chn(130,54-T_chn+t)+ep(130,t+1);
        end
        BB_chn(:,t+1)=BB_chn(:,t);
        PP_chn(:,:,t+1)=PP_chn(:,:,t);
        nm=find(M(:,54-T_chn+t)==1);
        DD=[diag(1-DD_chn(nm,54-T_chn+t)),diag(DD_chn(nm,54-T_chn+t))];
        PP_chn([nm;N+nm],[nm;N+nm],t+1)=PP_chn([nm;N+nm],[nm;N+nm],t+1)+DD'*(Sigma(nm,nm)\DD);
        Pd_chn=sqrt(diag(PP_chn(:,:,t+1)\eye(2*N)));
        BB_chn([nm;N+nm],t+1)=BB_chn([nm;N+nm],t+1)+PP_chn([nm;N+nm],[nm;N+nm],t+1)\(DD'*(Sigma(nm,nm)\(y_chn(nm,54-T_chn+t)-DD*BB_chn([nm;N+nm],t))));
        U0_chn=exp(alpha(130,1)+theta(1)*(Pd_chn(130,1)*sgi_nodes(:,1)'+BB_chn(130,t+1)+...
            ep_sd(130,1)*sgi_nodes(:,2)'));
        U0_chn=(U0_chn./(1+U0_chn))*sgi_weights;
        U1_chn=exp(alpha(130,1)+theta(2)*(Pd_chn(N+130,1)*sgi_nodes(:,1)'+BB_chn(N+130,t+1)+...
                ep_sd(130,1)*sgi_nodes(:,2)')-f(130,1));
        U1_chn=(U1_chn./(1+U1_chn))*sgi_weights;
        DD_chn(130,54-T_chn+t+1)=M(130,54-T_chn+t+1).*(U1_chn>=U0_chn);
    end
    D_chn(:,:,nsp)=DD_chn;
    B_chn(:,:,nsp)=BB_chn;
    P_chn=P_chn+PP_chn/NSP;
end
D_chn_ss=1*(mean(D_chn,3)>0.5);
B_chn_ss=mean(B_chn,3);

% for Figure 9
Pinv_chn=P_chn;
for t=1:size(P_chn,3)
    Pinv_chn(:,:,t)=P_chn(:,:,t)\eye(2*N);
end
unc_chn=zeros(N,size(P_chn,3));
for n=1:N
    for t=1:size(P_chn,3)
        unc_chn(n,t)=sqrt(Pinv_chn(N+n,N+n,t)+Pinv_chn(n,n,t)-2*Pinv_chn(N+n,n,t));
    end
end
writetable(table(unc_chn(130,:)),'../../Figures/Figure_9/china_unc_cfactual065.csv')


%% large shock (11%) to Mongolia, Taiwan, India, Japan, South Korea, Philippines, (2004-T_chn):2004
P_chn=zeros(size(P(:,:,54-T_chn:T-1)));
parfor nsp=1:NSP
    ep=epsilon(:,:,nsp);
    y_chn=y;
    y_chn([131:132,134:136,148],54-T_chn:54)=0.11+ep([131:132,134:136,148],1:T_chn+1);
    PP_chn=zeros(size(P(:,:,54-T_chn:T-1))); %#ok<*PFBNS>
    BB_chn=zeros(2*N,T-54+T_chn);
    % 2001
    nm=find(M(:,54-T_chn)==1);
    DD=[diag(1-D(nm,54-T_chn)),diag(D(nm,54-T_chn))];
    PP_chn(:,:,1)=P(:,:,54-T_chn);
    PP_chn([nm;N+nm],[nm;N+nm],1)=PP_chn([nm;N+nm],[nm;N+nm],1)+DD'*(Sigma(nm,nm)\DD);
    Pd_chn=sqrt(diag(PP_chn(:,:,1)\eye(2*N)));
    BB_chn(:,1)=B(:,54-T_chn);
    BB_chn([nm;N+nm],1)=BB_chn([nm;N+nm],1)+PP_chn([nm;N+nm],[nm;N+nm],1)\(DD'*(Sigma(nm,nm)\(y_chn(nm,54-T_chn)-DD*BB_chn([nm;N+nm],1))));
    U0_chn=exp(alpha(130,1)+theta(1)*(Pd_chn(130,1)*sgi_nodes(:,1)'+BB_chn(130,1)+...
        ep_sd(130,1)*sgi_nodes(:,2)'));
    U0_chn=(U0_chn./(1+U0_chn))*sgi_weights;
    U1_chn=exp(alpha(130,1)+theta(2)*(Pd_chn(N+130,1)*sgi_nodes(:,1)'+BB_chn(N+130,1)+...
            ep_sd(130,1)*sgi_nodes(:,2)')-f(130,1));
    U1_chn=(U1_chn./(1+U1_chn))*sgi_weights;
    DD_chn=D;
    DD_chn(130,54-T_chn+1)=M(130,54-T_chn+1).*(U1_chn>=U0_chn);
    for t=1:T-54+T_chn-1
        if DD_chn(130,54-T_chn+t)~=D(130,54-T_chn+t)
            y_chn(130,54-T_chn+t)=true_beta(130)*(1-DD_chn(130,54-T_chn+t))+true_beta(N+130)*DD_chn(130,54-T_chn+t)+ep(130,t+1);
        end
        BB_chn(:,t+1)=BB_chn(:,t);
        PP_chn(:,:,t+1)=PP_chn(:,:,t);
        nm=find(M(:,54-T_chn+t)==1);
        DD=[diag(1-DD_chn(nm,54-T_chn+t)),diag(DD_chn(nm,54-T_chn+t))];
        PP_chn([nm;N+nm],[nm;N+nm],t+1)=PP_chn([nm;N+nm],[nm;N+nm],t+1)+DD'*(Sigma(nm,nm)\DD);
        Pd_chn=sqrt(diag(PP_chn(:,:,t+1)\eye(2*N)));
        BB_chn([nm;N+nm],t+1)=BB_chn([nm;N+nm],t+1)+PP_chn([nm;N+nm],[nm;N+nm],t+1)\(DD'*(Sigma(nm,nm)\(y_chn(nm,54-T_chn+t)-DD*BB_chn([nm;N+nm],t))));
        U0_chn=exp(alpha(130,1)+theta(1)*(Pd_chn(130,1)*sgi_nodes(:,1)'+BB_chn(130,t+1)+...
            ep_sd(130,1)*sgi_nodes(:,2)'));
        U0_chn=(U0_chn./(1+U0_chn))*sgi_weights;
        U1_chn=exp(alpha(130,1)+theta(2)*(Pd_chn(N+130,1)*sgi_nodes(:,1)'+BB_chn(N+130,t+1)+...
                ep_sd(130,1)*sgi_nodes(:,2)')-f(130,1));
        U1_chn=(U1_chn./(1+U1_chn))*sgi_weights;
        DD_chn(130,54-T_chn+t+1)=M(130,54-T_chn+t+1).*(U1_chn>=U0_chn);
    end
    D_chn(:,:,nsp)=DD_chn;
    B_chn(:,:,nsp)=BB_chn;
    P_chn=P_chn+PP_chn/NSP;
end
D_chn_ls=1*(mean(D_chn,3)>0.5);
B_chn_ls=mean(B_chn,3);


%% for Figure 9

% mean beliefs
writetable(table(B(N+130,:)-B(130,:)),'../../Figures/Figure_9/china_original.csv')
writetable(table(B_chn_ss(N+130,:)-B_chn_ss(130,:)),'../../Figures/Figure_9/china_cfactual065.csv')
writetable(table(B_chn_ls(N+130,:)-B_chn_ls(130,:)),'../../Figures/Figure_9/china_cfactual11.csv')

% belief uncertainty
writetable(table(uncertainty(130,:)),'../../Figures/Figure_9/china_unc_original.csv')
Pinv_chn=P_chn;
for t=1:size(P_chn,3)
    Pinv_chn(:,:,t)=P_chn(:,:,t)\eye(2*N);
end
unc_chn=zeros(N,size(P_chn,3));
for n=1:N
    for t=1:size(P_chn,3)
        unc_chn(n,t)=sqrt(Pinv_chn(N+n,N+n,t)+Pinv_chn(n,n,t)-2*Pinv_chn(N+n,n,t));
    end
end
writetable(table(unc_chn(130,:)),'../../Figures/Figure_9/china_unc_cfactual11.csv')



